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ABSTRACT 

Wc present a new pseudospectral algorithm for the calculation of the structure of 
atoms in strong magnetic fields. Wc have verified this technique for one, two and three- 
electron atoms in zero magnetic fields against laboratory results and find typically 
better than one-percent accuracy. We further verify this technique against the state- 
of-the-art calculations of hydrogen, helium and lithium in strong magnetic fields (up 
to about 2 X 10^ T) and find a similar level of agreement. The key enabling advantages 
of the algorithm are its simplicity (about 130 lines of commented code) and its speed 
(about 10^ — 10^ times faster than finite-element methods to achieve similar accuracy). 



Key words: physical data and processes: atomic data — physical data and processes: 
magnetic fields — stars: neutron — stars: magnetic fields — stars: white dwarfs 



1 INTRODUCTION 

A hydrogen atom in a magnetic field is arguably the sim- 
plest Hamiltonian without an analytic solution - it is a com- 
bination of a harmonic oscillator with a Coulomb potential. 
The problem of atoms in magnetic fields is also a vexing 

■ challenge in the study of neutron stars and white dwarfs. A 

■ sufficiently intense magnetic field cannot be treated pertur- 
, batively. If the magnetic field completely dominates over the 

■ nucleus one can un derstand the a tomic structure as a one - 
' dimensional atom (|Loudonlll959l : [Haines fc Robert^ 1 1969h . 
, A non-trivial relaxation of the one-dimensional assumption 

■ is the adiabatic approximation that assumes that the wave- 
function of the electron perpendicular to the magnetic field 
is simply the gr ound Landau level of th e electron without 
the nucleus (e.g. iHevl fc Hernauist|[l998l ). 

The most challenging regime is of course where the nu- 
cleus and the magnetic field compete for the electron's atten- 
tion. This lies around 10'*^" T for hydrogen and is known as 
the strong-field regime. For such field strengths one cannot 
assume that the electron lies in the lowest Landau level. For 
a detailed bibliography of the vast literature of atomic cal- 
culations in strong magneti c fields, we enco urage the reader 
to consult our recent work ([Thirumalai fc H cvl 20C)!|). 

In the strong-field regime, one can assume that the 
wavefunctio n is a sum of La ndau levels or spherical har- 
monics (e.g. iRuder et al.lll99^ ) or solve the two-dimensional 
partial differential equat ion without this assumption (e.g. 
iThirumalai fc Hevll 12003 '). This work takes the second ap- 
proach and solves the differential equations for the electronic 
wavefunctions and the interelectron potentials numerically. 



The main goal of this paper is not to present more accurate 
calculations of atoms in strong magnetic fields (although the 
calculations rival the current state of the art), but rather the 
purpose is to provide a new, much faster algorithm to per- 
form these calculations and verify its results against previous 
work. 

This new algorithm is straightforward to implement and 
compact; the thirty-line program for single-electron atoms 
and the 130-line program for multiple-electron atoms are 
presented in the appendices. With this simplicity comes a 
dramatic speed-up of the calculation o f a factor of 10^ — 
10^ r elative to finite-element techniques (|Thirumalai fc HevJ 
l2009h while achieving similar accuracy. This pseudospectral 
algorithm interpolates the approximate solution to the dif- 
ferential equations as a polynomial over an irregular mesh. 
Because the approximate solution is a polynomial, it is also 
analytic over the mesh. Even if the actual solution oscillates 
on a scale of a couple of mesh points, the analytic interpolant 
still provides an accurate approximation to the solution. On 
the other hand, if the actual solution develops discontinu- 
ities (i.e. it is not analyt ic), the pse udospectral method is 
less useful, (however see lBovdll200ll ). In the case of mag- 
netized atoms, the solution to the Schrodinger equation is 
analytic everywhere except at the origin where we can set 
a boundary condition, so the pseudospectral method yields 
great returns. 

In the next sections we will present the partial differen- 
tial equations for the electronic wavefunctions ( i]2.ip . the in- 
terelectron potential ( tj2.2p . an introduction to pseudospec- 
tral methods f i^2.3p . the application of boundary conditions 
( i]2.4l) . our results (^ and source code (Appendices RllB)) . 
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2 CALCULATIONS 
2.1 The Hamiltonian 

Rather than examine the Hamiltonian of an atom in a mag- 
netic field in detail, we will simp ly recap the relevant equa - 
tions here and refer the reader to lThirumalai fc Hevll (|2009l ) 
for further details. We use the Hartree-Fock equations to de- 
scribe the wavefunction of a multiple-equation system, this 
yields a coupled set of Schrodinger-like equations for an ap- 
proximation to the wavefunction of each electron. 



rs Z 



(1) 



where i, j — 1, 2, 3, 



,iV. 



Z is the charge of the nucleus located at the origin, nj 



Here A'^ is the number of electrons, 

Vij is 

the inter-electron potential operator and Pz measures the 
strength of the magnetic field, /3z = Beh / {2Z^ m^c?) = 
B/(4.7Z^ X lO^'T). The variable is the distance from the 
origin and R± = Ts smO, the distance from the ax;is defined 
by the direction of the magnetic field; this notation is con- 
sistent with the software in the appendices. The functions 
must be finite at the origin and approach zero as r — >■ oo. 
The total binding energy of the atom is given by 



E ■ 



E 



E, 



rrii + Si 



E 



{■4}i\Vij\tpi 



(2) 



in units of Z^a^rrisf? /2 

Although the Hamiltonian is no longer spherically sym- 
metric, it is still symmetric about the magnetic axis, so we 
expect the azimuthal or magnetic quantum number rrii to 
be good. We will label the solutions to Eq. ^ by Vi, Si and 
mi where Si is the spin of the electron and Ui equals 1 for 
the most bound state with a particular value of Si and rrii 
and 2 for the next most bound state and so on. 

The first step is to specify the precise differential equa- 
tions to solve. Specifically, we will take 



tpi 



(3) 



The parameter ^ = 008 6^, where 6 is the angle from the z- 
axis. For the radial dependence, the domain is [0, co), so we 
would like compactify the domain to [—1, 1]. Let 



2 tanh ■ 



rz 



(4) 



where ranges from zero to infinity while r now ranges 
from —1 to 1. The quantity rz is the "zoom radius" that 
determines at what value of rs the relationship between r 
and Ts goes from linear to exponential. 



2 1, r <C rz; r ! 

rz 



1 - 4exp(-2rs/rz),r > rz. (5) 



Because 2 tanh 1 — 1 pa 0.5 about one third of the points 
lie at r > rz. Because the domain is compactified we can 
easily apply boundary conditions both at infinity and at the 
origin as appropriate for solving the Schrodinger equation: 
iti( — = Ui(-fl,/i) = and the additional condition for 
mi / that Ui{r, — 1) = Ui{r, +1) = 0. 

Using the substitution of Eq. ((3| into ([1]), yields 



9ri 



I s 



r| 9/i 



imi4i 



where 
(9r| 

where r^ = (r -I- l)/2. 



4(1--^) 



u, = 



(6) 



(7) 



2.2 Inter-electron potential 



Using the results of H2.3I we can convert this partial differ- 
ential equation into a matrix equation, specifically an eigen- 
value equation. However, before doing this in H2.41 we must 
discuss the inter-electron potential Vij. 

The inter-electron potential consists of two terms, the 
direct and exchange interactions where 



Vdirc 



Vc 



exchange, 



(8) 



if the spins of electrons i and j are aligned and Vij = Vdircct.ij 
otherwise. The direct interaction satisfies 



Vd\vcct,ij1pi = (j}j{v)'4)i 

where 



-4:TTtp*tpj 



(9) 



(10) 



which we can solve, knowing the values of xpj, by inverting 
the Laplacian to yield 



0. = (V) ' [-4^i/.*^i 



(11) 



where (V^) is the inverse of the Laplacian with the appro- 
priate boundary conditions supplied. Therefore, the direct 
interaction is a diagonal operator, essentially a function. 

On the other hand the exchange interaction is more 
complicated. We have 



^exchange , ij i^i = 4>ij {r)i}j 



where 



-4TT^p*ll)i 



(12) 



(13) 



Therefore, we have 

Voxchangcij'/'i = (^ij(r)?/)j = (V^) ^ [-471?/)* l/^i] . (14) 

Thinking of the different components of the right-hand sides 
of Eqs. (|12p and (|15|) as matrices, we can write the poten- 
tials due to the direct and exchange interactions in matrix 
notation as 



Kxchangci = -47r diag(?/>j) (V^) diag(V'j) 



(15) 



where we have dropped the index i because the exchange 
interaction operator only depends on the wavefunction of 
the electron j. Both the direct and exchange potentials must 
go to zero at r ^ oo and be finite at the origin. 

The set of coupled equations is of course non-linear and 
difficult to solve directly; however, we can approach the so- 
lution iteratively, by first solving the equations ([1} ignoring 
the inter-electron terms and then on subsequent iterations 
using the wavefunctions from the previous iteration to cal- 
culate the interaction operator using Eq. (Illf) and (|15|l . At 
each iteration the equations are decoupled, so the time to 
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calculate a single iteration scales as the number of electrons, 
making the calculation of even large atoms tractable. If the 
equations were coupled through the exchange term, then 
each dimension of the matrix to diagonalize would scale with 
the number of electrons so the time to complete an iteration 
would scale as the cube of the number of electrons. 

Let's evaluate for the exchange potential operator using 
the substitution of Eq. Q. We have 



exchange, 



: <i)ij(r)uj 



and 



• exchange, ij 



Ui = (t>ij{v)uje 



(16) 



(17) 



where Am = rm — mj . Because none of the other terms in 
the Hamiltonian depend of neither side of Eq. (|17p can 
depend on (j) so let us define 



;(r) = ij)ij"'{rs,^i)t 



where 



Simplifying this result yields 



'(rs,lJi) = -47r- 



(18) 
(19) 

(20) 



and 



-47r diag(itj) ( V 



,2 Am^ 



diae 



where 



+ 



1-M 



(21) 



(22) 



The result for the direct interaction is nearly trivial be- 
cause it is a function rather than an operator. We have 



rectjij 



= -4^(V^)"' 



(23) 



We must specify the boundary conditions on the inverse 
operators in Eq. (|21l) and 1)23^ . Both the direct and exchange 
potentials must go to zero at r — >■ cxd and be finite at the 
origin. Furthermore, for Am 7^ 0, the exchange potential 
must vanish along the magnetic axis. 

In principle one could solve for all of the electronic wave- 
functions at the same time, but this would render the dif- 
ferential equations non-linear; therefore, we first neglect the 
electronic contribution to the potential to find a first guess at 
the electronic wavefunctions. These wavefunctions become 
the input into Eq. \21\ and (|23p for the second iteration. 
This process is repeated until the total energy of the state 
is constant to one part in 10^ between successive iterations. 



2.3 Pseudospectral methods 

The basic idea behind interpolating spectral (or pseudospec- 
tral) methods is to approximate functions by their values at 



a set of points and a set of smooth functions to interpolate 
between these points for the purpose of calculating deriva- 
tives and integrals. An example is useful to make this con- 
crete. A function f{x) takes the values fi at the points Xi. We 
could construct an interpolation scheme by linearly interpo- 
lating between the points or using a cubic spline. This is the 
basis of finite difference schemes with successively higher ac- 
curacy. However, if we wanted an analytic function through 
the points we could use the polynomial that passes through 
all of the points. This polynomial interpolant is given by 



(24) 



where we have used the notation of iTrefethenI (|2000l ). The 
polynomial Pi{x) has the property that it takes the value of 
unity at Xi and it vanishes at all of the other points where 
the function values are given. We can take the derivative of 
the basic polynomial to yield 



Pi{x) =pi{x)^{x - Xj) ^ 

and evaluate it at one of the points a:*; to yield 
Pi{xk) = Pi{xk) ^(a;fe - Xj)~'^ = Afe 

where 



Afe = t4 / = y ^(a 



(25) 



(26) 



(27) 



where Dik is the differentiation matrix such that ff{xk) ~ 
Dikfi- Unlike finite element, finite volume and finite differ- 
ence techniques, the differentiation matrix is dense, so it is 
more costly to manipulate; however, one can often achieve 
similar accuracy with much smaller matrices for problems 
with smooth solutions. This decrease in the size of the ma- 
trices more than offsets the burden of handling dense matri- 
ces. 

The choice of the points Xi is completely arbitrary but 
for finite domains it is more stable to pick the points to be 
bunched to ward the ends o f the domain and more sparse in 
the middle (|Trefethenll2000l ) . A convenient set of points over 
the domain [—1, 1] are the Chebyshev points 



cos{iTv/N),i = 0...N. 



(28) 



This domain is well suited for the angular dependence in 
spherical coordinates where we define /i = cos^ — z/r and 
use the Chebyshev points to sample /i. Thus the angular 
dependence is evenly sampled in 6. 



2.4 Discrete equations 

The Octave programs given in the appendices for the single 
and multiple-electron problems convert the partial differ- 
ential equations of i)2.1l and 12.21 to matrix equations using 
the pseudospectral methods outlined in t |2.3l We encourage 
readers that are familiar with Octave or Matlab to look un- 
der the hood and dissect the routines. 

The previous section outlined how to construct a differ- 
entiation matrix using a pseudospectral method in one di- 
mension. The second derivative is simply the square of this 
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matrix and an integration operator is the inverse of the dif- 
ferentiation matrix with some boundary condition supplied. 
The remaining subtleties are how to apply boundary con- 
ditions and how to move to more than on e dimension. The 
techniques outlined here are presented in iTrefethenI (|200(]| ') 
in greater detail, clarity and generality. 

Supplying the boundary conditions for the eigenvalue 
equations is rather straightforward and most easily illus- 
trated by an example. Let us find the eigenvalues of the 
matrix M with the solution / vanishing at the end points. 
We have 



Mf = A/. 

Let's define the diagonal matrix, 



(29) 



B = 




1 









(30) 



with all of the diagonal elements unity except for the first 
and last that vanish. Multiplying the vector / by _B has the 
effect of setting the endpoints to zero, so let's look at the 
following eigenvalue equation 



BMBf = A/. 



(31) 



The matrix B that precedes M ensures that the resulting 
output vector satisfies the boundary conditions, and the ma- 
trix that follows M ensures that the input vector satisfies 
the boundary conditions as well, so the eigenvectors of the 
matrix BMB except for two are guaranteed to satisfy the 
boundary conditions. The two remaining eigenvectors are 
linear combinations of [1 ■ ■ ■ 0] and [0 ■ • ■ 01] and have 
zero eigenvalues. We can either ignore these eigenvectors or 
define a new matrix M that lacks the first and last columns 
and the first and last rows of the original matrix M. The 
eigenvectors of Af are the same at those of BMB but with 
the first and last entries omitted. 

In summary, to obtain the set of eigenvectors of M that 
satisfy the boundary conditions we construct M and add 
the first and last entries if needed. For the potential calcu- 
lation it is often more useful to work with BMB + {I — B) 
so that the resulting matrix can be inverted before setting 
the boundary conditions. In the potential calculation we set 
the boundary condition at infinity because the Laplacian 
matrix vanishes for the row and column corresponding to 
infinity, so we drop the row and column before inverting the 
Laplacian. However, we would also like to set a boundary 
condition at the origin, so we must include the origin in the 
calculation. At the origin the angular terms in the Laplacian 
diverge, but this is a coordinate singularity that we address 
by including only a single entry for the origin and drop the 
angular derivatives there. In this way, we address the diver- 
gent terms at the origin and the vanishing terms at infinity. 

The second issue is how to build the various operators 
in more than one dimension. Let's say that we have obtained 
a vector of points x of dimension and a pseudospectral 
derivative Dx and similarly a vector of points y of dimension 



where Ijv is a vector of A'^ ones and (8> denotes the Kronecker 
product, we obtain a two-dimensional mesh {x' , y) of points 
in two vectors of dimension N x M. Finally we would like 
to construct the Laplacian over this mesh (this example is 
Cartesian). We have 

^2 



'lM+lN®{DyY 



(33) 



1m, y = In ®y 



(32) 



where Im is the N x A'-identity matrix. The resulting matrix 
has dimensions (A' x M) x (A x M) and calculates the 
Laplacian over the mesh specified by [x ,y'). 

These techniques allow us to discretize the equations 
outlined in tj2.ll and l2.2l to obtain the software in the appen- 
dices. 



3 RESULTS 
3.1 Hydrogen 

As a first test of the algorithm, we examined hydrogen in 
a strong magnetic field and in zero field. The Octave pro- 
gram is presented in Appendix [X] This is one case where 
the program is shorter than the tables it calculates and re- 
places. In zero-field with 30 radial points (A = 31, exclud- 
ing excluding —1 and 1 in the compactified domain, corre- 
sponding to the origin and infinity, respectively in the real 
domain) and with 12 angular points (M = 11), it calculates 
the eigenvalues up to principal quantum number, n = 12 
to at least three digits of accuracy (a total of 650 states). 
Up to /3 ~ 10 it determines the low-lying states of hydrogen 
(again with A = 31 and M = 11) to at least three digit s 
of accuracy (|Thirumalai fc Hevll [20091 : iRuder et aLlll994l '): 
for /3 1, the accuracy is typically two orders of magni- 
tude better. The distinct advantage of the pseudospectral 
algorithm is the spee d to achieve accurate resu lts. The fine 
mesh calculations of I Thirumalai fc Hevll (|2009h took up to 
two days or 1.7 x 10^ seconds — the pseudospectral algo- 
rithm takes about 1.4 s to get the same accuracy on the 
same processors. The algorithm uses spherical coordinates 
so it is not well suited to fields much stronger than /3 = 10. 
At /3 = 100, a finer mesh (M = 31 and A = 51) and 150 
seconds of computation t ime are requ i red to obtain three 
digits of agreement with I Ruder et al.l (| 19941 ): however, in 
this regime, cylindrical coordinates may be more appropri- 
ate and one can use the adiabatic appro ximation to obtain 
accurate results (|Hevl fc Hernauis3ll998l ). 



3.2 Helium and Lithium 

The problem of multi-electron atoms (Appendix |B]) is of 
course much more complicated than hydrogen. Again we 
begin with the zero-field case, but we have many states to 
examine that test the direct and exchange interactions sep- 
arately. In this case some validation for a vanishing mag- 
netic field is useful. Table. [1] compares the zero-field re- 
sults for some low-lying states of helium and lithium with 
the experimentally determined values from iRalchenko et al.l 
(2008). Again the algorithm achieves typically better than 
one-percent accuracy rapidly. 

Of course, the technique is not intended to supplant the 
p opular, a c curat e and fast MCHF algorithms such as those 
of lFischeil j 19971 ). Such algorithms achieve their speed and 
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Table 1. Compa rison of the observed z ero-field states of helium 
and lithium from lRalchenko et al ] l|2008l 'l to the values calculated 
here. Energies arc in units of Rydberg energies in the Coulomb 
potential of nuclear charge Ze, where Z = 2 for helium and Z = 3 
for lithium. We use M = 7, = 31 and rz = 58. The value of rz 
is chosen to provide close agreement for the ground state of helium 
and held fixed. Individual calculations take about ten seconds for 
helium and twenty for lithium. The primary sourc e for the helium 
data is lDrake fc Martini l l 19981) and lMoorj l ll97ll ) for the lithium 
data. 



Conf. 


Term 


J 








A(%) 


Helium I 









1 


451692958 


1.43189 


-1.364 


ls2s 




1 


1 


087514256 


1.08996 


0.225 


ls2s 







1 


072885081 


1.07929 


0.598 


ls2p 




2 


1 


066484964 


1.06856 


0.195 






1 


1 


066484790 


1.06857 


0.196 









1 


066482539 


1.06856 


0.195 


ls2p 


Ipo 


1 


1 


061818980 


1.06593 


0.387 






1 


1 


061818980 


1.06594 


0.388 


ls3s 




1 


1 


034248834 


1.03719 


0.284 


ls3s 







1 


030539891 


1.03456 


0.391 


ls3p 


3po 


2 


1 


028945784 


1.03174 


0.272 






1 


1 


028945735 


1.03174 


0.271 









1 


028945117 


1.03174 


0.272 


ls3d 




3 


1 


027722445 


1.03073 


0.293 






2 


1 


027722439 


1.03073 


0.293 






1 


1 


027722338 


1.03073 


0.293 


ls3d 


ID 


2 


1 


027714652 


1.03073 


0.294 


ls3p 




1 


1 


027476815 


1.03092 


0.335 


Lithium I 


ls^2s 




1/2 


1.661773620 


1.65595 


-0.350 


ls^2p 


2po 


1/2 


1.646683387 


1.64097 


-0.347 






3/2 


1.646683042 


1.64097 


-0.347 


ls^3s 


25 


1/2 


1.634226909 


1.62883 


-0.330 



accuracy by using the spherical symmetry of the atom: for 
example, by using spherical harmonics and a special treat- 
ment of full shells. The power of this pseudospectral algo- 
rithm becomes crucial when spherical symmetry is destroyed 
by a strong magnetic field. 

Table. [2] traces the binding energy of two low-lying 
states of helium using the soft ware in Appendix [B] The 
agreement with our earlier work (|Thirumalai fc Hevlll200^ ) 
is striking at weaker fields but as the magnetic field starts to 
dominate the atom, the agreement is poorer. Better agree- 
ment is achieved by increasing the angular resolution of the 
calculation to account for the changing shape of the atom. 

For lithium we comp are with the results of 
lAl-Huiai fc Schmelcheij \2004 ) who use a different pa- 
rameter for the strength of the field, 7 — 2Z^ /3z = 18/3z 
where the second equality holds for lithium. We have chosen 
to use eight angular grid points and thirty-two in the radial 
direction. 

Here the agreement with the previous work is encour- 
aging at least up to moderate field strengths (7 < 1). For 
stronger fields the pseudospectral method typically yields 
excessively bound systems, except for the 3)+ state 
that happens to be the most bound state in strong mag- 
netic fields. These spurious eige nvalues are not uncommon 
for pseudospectral methods (e.g. lBovdll200ll ). The software 
as presented in the appendices actually includes a filter for 
a blatantly spurious eigenfunction that attains a large value 



Table 2. Table listing the binding energies of the most tightly 
bound states of helium in moderate to large magnetic fields; = 
— 1 or M2 = —2 and Sz = —1,iTz = +1. Energies are in units of 
Rydberg energies in the Coulomb potential of nuclear charge Ze, 
where Z = 2 for helium. The results fro m the current work can be 
compared readily with previous work llThirumalai fc Hev]||2009l ; 
iRuder et aT]|l994l) . The first row for /3z = 1 is for M = \\ like 
the other rows; the second row gives the result s from M = 41. 
TH08 indicates result s fromlThiruraalai fc Hevll ||2009|) . and R94 
indicates results from lRuder et ah I ||1994| '|. 







Mz = -1 






Mz = -2 






Here 


TH08 


R94 


Here 


TH08 


R94 


0.01 


1.1221 


1.1183 


1.1182 


1.0871 


1.0852 


1.0828 


0.02 


1.1656 


1.1612 


1.1609 


1.1265 


1.1234 


1.1207 


0.05 


1.2734 


1.2691 


1.2658 


1.2215 


1.2175 


1.2097 


0.07 


1.3359 


1.3319 


1.3258 


1.2764 


1.2732 


1.2596 


0.10 


1.4211 


1.4189 


1.4069 


1.3515 


1.3510 


1.3266 


0.20 


1.6586 


1.6585 


1.6270 


1.5653 


1.5598 


1.5073 


0.50 


2.1603 


2.1550 


2.0508 


2.0415 


2.0009 


1.8508 


0.70 


2.4584 


2.4029 


2.2329 


2.2902 


2.2246 


1.9935 


1.00 


3.2165 


2.7026 


2.4675 


3.0102 


2.4981 


2.2572 


1.00 


2.7115 






2.6668 







Table 3. Table listing the binding energies of three tightly 
bound states of lithium in moderate to large magnetic fields. 
At the top of each pair of columns is the symmetry subspace 
2S+l(;;^^-)z-parity -y^g ^^y-^ ^^le dominant configuration for the 
ground state wi thin each subspace. Th e column "Other" gives 
the values from iRalch enko et al\ l|2008l ) for zero field and from 
lAl-Huiai fc SchmelcheJ l l2004h elsewhere. 



7 


•2 


0+ 


H- 


1)+ 




-3)+ 


18/3z 


Here 


Other 


Here 


Other 


Here 


Other 


0.00 


1.6560 


1.6618 


1.1968 


1.1925 


1.1358 


1.1299 


0.01 


1.6590 


1.6629 


1.2017 


1.1969 


1.1427 


1.1487 


0.05 


1.6681 


1.6673 


1.2192 




1.1655 


1.1663 


0.10 


1.6775 


1.6705 


1.2390 


1.2334 


1.1901 


1.1869 


0.20 


1.6928 


1.6741 


1.2736 


1.2674 


1.2330 


1.2278 


0.50 


1.7255 


1.6729 


1.3533 


1.3463 


1.3362 


1.3294 


1.00 


1.8256 


1.6575 


1.4821 


1.4432 


1.4713 


1.4627 



along the magnetic axis for large distances from the nucleus. 
It is our experience with the hydrogen atom that these spu- 
rious eigenvalues become more prevalent for stronger mag- 
netic fields and are more bound than the eigenvalues for well 
behaved eigenfunctions. However, the algorithm still finds 
the correct eigenvalue for the ground state, so with a more 
sophisticated filtering technique spurious eigenvalues can be 
eliminated for calculations in higher field strengths than con- 
sidered here. 



4 DISCUSSION 

One of the primary goals of the current work was to pro- 
vide a computational method that would be rather econom- 
ical with regard to computation time, without having to 
compromise on accuracy. As can be seen in Table [1] the 
desired level of accuracy can be achieved with computa- 
tion times on the order of seconds. The calculated values 
therein can conceivably have better agreement with experi- 
mentally determined values if effects of spin-orbit coupling, 
relativistic corrections and the effects of mixing of config- 
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urations were to be included. Similarly, upon examining 
Table [2j it is immediately apparent that the estimates of 
the binding energies of the state of helium corresponding to 
Mz = —l\Sz = —1 are more bound than the previous esti- 
mates. The average improveme nt with respect to the val ues 
reported in our previous study (jThirumalai fc Hevlll2009l ). is 
about 0.5% while, the average improvement over the values 
reported in lRuder et all j 19941 ) is about 3.4% over the entire 
range of magnetic field strengths reported therein. Similarly, 
upon comparing the results from the present calculation for 
the state of helium corresponding to = —2;Sz = — 1 
we see that they are more bound on average by about 1.5% 
and 5.8% respectively, w ith r espect to the calcu lations of 
iThirumalai fc Hevll |200i) and lRuder et all (|l994h . Finally, 
when comparing with the results of other authors in Table [3] 
we see that the current single-configuration calculation us- 
ing spectral methods provides binding energies that are an 
improvement on previous estimates for three tightly bound 
states of lithium. It can be seen that the average range of 
improvements for the estimates of the binding energies of 
the three states, relative to the calculations of other work, 
is about 1%. 

In carrying out the calculations in the current study, 
as was noted earlier, the domain was compactified accord- 
ing to Eq. Q. This has the distinct advantage that one 
can set a boundary condition at true infinity rather than 
at a large value of the radial distance r. Typically, such a 
procedure renders the partial differential equations highly 
non-linear, making their solution with finite-element tech- 
niques more involved and computationally demanding, thus 
requiring a greater amount of computation time. In the cur- 
rent study, the compactification effectively changes the op- 
erators in Eqs. (O, ([21} and pSjl . 

Since the algorithm presented here is essentially a sim- 
ple prototype, several avenues for further development are 
apparent. They address the stability, accuracy and robust- 
ness of the calculation. First from a physical point of view it 
does not make sense to use the same zoom radius (the value 
of rz) for all of the electrons in an atom. The zoom radius 
is essentially the scale within which the radial mesh is fine; 
therefore, it should reflect the expected or calculated extent 
of the electronic wavefunctions. It would be quite natural 
to adaptively change the value of rz to be some constant 
factor times (r) or (r^) as the Hartree-Fock iteration pro- 
ceeds. Because the interelectron potential is smoother than 
the underlying wavefunctions it would be more effective to 
calculate the potentials for each electron using its own value 
of rz and then interpolate these results onto the meshes of 
the other electrons. In this manner the value oi rz would 
adapt for the particular orbital that the electron occupies, 
achieving higher accuracy without increasing the number of 
radial mesh points. For the direct interaction this is straight- 
forward because it is essentially a function; however, the ex- 
change interaction is an operator or a matrix, as in Eq. psp . 
so the interpolation of the final result is a bit problematic. 
Rather it is expedient to interpolate the wavefunctions of 
the other electrons onto each other's mesh and calculate the 
exchange interaction one pair at a time. Although this in- 
terpolation would necessarily make the code more cumber- 
some, it need not increase the processing time significantly. 
The inverse Laplacians for different values of rz are related 
to each other through simple scalings so the costly matrix 



inversions would still only happen once at the beginning, 
and the results would scale for the appropriate values of rz- 
A similar adaptation can be performed in the angular 
direction to focus grid points along the magnetic axis for 
strong magnetic fields. In this way we can achieve high ac- 
curacy for all of the electrons without including more grid 
points. A simple way to achieve this is the "Arctan/Tan" 
mappmg (|Bovdll200ll ') where the original angular variable 6 
is mapped onto 



usmg 



arctan (L tan ^ 



(34) 



The parameter L can change to refiect the actual physical 
extent of the wavefunction in the angular direction for each 
electron in the atom as the iterative process proceeds. This 
adaptation would add little additional overhead above the 
radial adaptive scheme outline in the last paragraph. 

A second issue that we encountered is the appearance 
of spurious eigenvalues among the physical ones, especially 
for intense magnetic fields. The current software eliminates 
these eigenvalues by looking at the structure of the eigen- 
functions. Essentially, many of the spurious eigenfunctions 
have significant support along the magnetic axis at large 
radii. Unfortunately, the current set of criteria are not suffi- 
ciently robust to exclude all of the spurious eigenfunctions. 
There are two clear paths to address this issue. The first is 
to develop a more sophisticated filtering technique to check 
the eigenfunctions during the iterative process. The second 
technique would be to use cylindrical coordinates that are 
better adapted to the intense-field limit where the wave- 
functions essentially are well approximated by the product 
of Landau orbitals and a function of the coordinate along 
the magnetic field (the so-called adiabatic approximation). 
The use of cylindrical coordinates is of course not equiva- 
lent to the adiabatic approximation but it does fit better 
with the dominant symmetry of the intense-field problem. 
For this work we did not employ cylindrical coordinates be- 
cause we wanted to verify our results against the zero-field 
limit where the binding energies of the various species are 
known precisely. 

From a mathematical point of view, we have chosen to 
use Chebyshev points to sample in both the angular and 
radial directions. In this particular case it is possible to ex- 
ress the spectral derivative s using a fast-Fourier transform 
Trefethenllioool : lBovdll2"00ll ) . The most obvious application 
for the problem at hand would be to speed up the calcu- 
lation of the direct component of the inter-electron poten- 
tial; because this only makes a modest contribution to the 
workload, we decided to forgo the added complexity and use 
matrix multiplication to determine the inter-electron poten- 
tial. On the other hand, given some additional work it may 
be possible to speed up other parts of the calculation using 
FFTs. 

A final improvement would be to include sev- 
eral electronic configuratio ns for each state (e.g. 
lAl-Huiai fc Schmelchei) |2004| ) not just the dominant 
one as we did here. Given the excellent agreement with the 
observed atomic properties that we found for the zero-field 
limit, we expect the improvement of a multi-configuration 
calculation to be modest at least for small atoms. However, 
recent theoretical results indicate that the atmospheres of 
neutron stars may be composed of elements near silicon 
(IHoffman fc Hevll |2009| ). so including many configurations 
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may be crucial at least to verify the algorithm in the 
weak-field limit. 



5 CONCLUSION 

The amount of literature available regarding atoms in strong 
magnetic fields is quite extensive. Since the 1970s with im- 
provements in both computing infrastructure and numerical 
methods the problem of atoms in strong magnetic fields has 
become increasingly tractable. However, the current state- 
of-the-art calculations still take on the order of hours if not 
days to obtain results of reasonable accuracy. To address 
this situation, we have, in this study, presented a new pseu- 
dospectral method for the calculation of the structure of 
atoms in strong magnetic fields. The key enabling advan- 
tages of the algorithm are its simplicity (about 130 lines 
of commented code) and its speed (about 10^ — 10^ times 
faster than finite-element methods to achieve similar accu- 
racy). For hydrogen, helium and lithium it gives results that 
agree with the previous work at the few percent level or bet- 
ter for fields weaker than about 2 x 10*" T. 

For the purpose of analyzing the spectra of magne- 
tized white dwarfs and neutron stars it becomes crucial 
to have accurate data for the energy levels of different 
atoms in strong magnetic fields. With the ability to per- 
form atomic structure calculations with a significant reduc- 
tion in computation time, as with the pseudospectral meth- 
ods described in this study, it becomes possible to amal- 
gamate atomic structure calculation software with atmo- 
sphere models of neutron stars and white dwarf stars. This 
is a direct advantage of the software provided in the ap- 
pendices to this paper. We have additionally presented sev- 
eral avenues for further research with the software includ- 
ing larger atoms, multi-configuration Hartree-Fock and full 
configuration-interaction calculations. 
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APPENDIX A: MAGNETIZED HYDROGEN 

Both the single-electron and multiple-electron codes are self-contain ed Octave prog rams with the exception of 
the eight-line function cheb that Trefethen published previously in IXrefethenl ()200C ) and is available online at 
|http : //www. comlab. ox . ac.uk/nick. tref ethen/cheb .m This function calculates the differentiation matrix in Eq. H27p . Both 
programs are also compatible with Matlab with some minor modifications in the output statements in the multi-electron 
program. However, after experimentation we concluded that the Octave version ran faster with our setup. 

The hydrogen code is written as a function to be called from the Octave command line that returns the eigenfunctions, 
eigenvalues, angular and radial points, a list of the eigenfunctions that pass a rudimentary test and the "zoom" factor. For 
values of M and A'^ on the order of ten, the vast majority of the computation time is devoted to calculating the eigenvalues, 
so the total execution time depends on the efficiency of the eigenvalue routine (eig); consequently, porting this algorithm to 
C does not result in a substantial speed up — in fact in our experience. Octave runs faster. 

7, HYDROGEN_S - calculate hydrogen in spherical coordinates 
7, beta is B/Bc, M is number of mu, N is number of radial 

the wavefunction is phi=f (r, mu=cos (theta) ) /r exp(i m phi) 
7. eigenfunctions in V, e-val in Lam, coords in (w,rs) 
function [V, Lam, w, rs , igood, zoom] = hydrogen_s(beta,M,N,mphi) 
7. w coordinate, ranging from -1 to 1 

[Dw,w] = cheb(M); D2w = Dw"2; 

if (mphi~=0) w=w(2:M); Dw=Dw(2 :M, 2 :M) ; D2w=D2w(2:M,2:M) ; end 

Hw=-diag(l-w. *w) *D2w+diag(2 . *w) *Dw; 
7. r coordinate, ranging from 1 to -1, rp from 1 to 

[D,r] = cheb(N); rp = 0.5*(r+l); D = 2*D; D2 = Y)'2; 

hh = diag(l-rp.*rp) ; D2 = hh* (hh*D2+diag(-2*r) *D) ; D = hh*D; 

D2 = D2(2:N,2:N) ; D= D(2:N,2:N); rp=rp(2:N); 
7. zoom factor: set by coulomb and larmor radius; , rs from inf to 

zoom=l/(l/110+sqrt(beta)/41) ; rs=zoom*atanh(rp) ; R = diag(l./rs); R2 = diag(l . /(rs . *rs) ) ; 

Hr=-l/(zoom*zoom)*D2-2*R; [rr.ww] = meshgrid(rs,w) ; rr = rr(:); ww = ww(:); rperp2=rr . *rr . * (1-ww. *ww) ; 

if (mphi==0) 

H = kron(Hr,eye(M+l))+kron(R2,Hw)+diag(beta*beta*rperp2) ; 
else 

H = kron(Hr , eye(M-l) )+kron(R2 ,Hw)+diag(beta*beta*rperp2) +diag(mphi*mphi . /rperp2) ; 
end 

[V,Lam] = eig(H) ; Lam = diag(Lam) ; [Lam,ii] = sort(Lam); Lam = Lam+2*beta* (mphi-1) ; V = V(:,ii); 
7i check outer B.C. and for bound states 

igood = find((V(l, :) .*V(1, :)) '<(M*N)-(-2)*le-4 & Lam<0) ; 



APPENDIX B: MAGNETIZED ATOMS 

The atomic code is written as a script to be called from the Octave command line because it has many parameters and 
diagnostics that one might like to examine and the script format encourages the user to choose what he or she wants to use. 

Similar to the case of hydrogen, for values of M and on the order of ten, the vast majority of the computation time 
is devoted to calculating the eigenvalues and inverting the Laplacians, so the total execution time depends on the efficiency 
of the eigenvalue and matrix inversion routines (eig and inv); consequently, we do not expect that porting this algorithm 
to a lower-level language [e.g. C) will result in a substantial speed up. In fact the computational effort is dominated by the 
diagonalization that is performed for each electron at each iteration. The inversion of the Laplacian is performed before the 
iterations begin. 

7. THREEDATOM - calculate an atom in three-D 
7. the wavefunction is phi=f (r,mu=cos (theta) ) /r 
7« r coordinate, ranging from -1 to 1 
7o Examples for Helium I 



nu= 


[1,1] 


spin= 


[+1/2,-1/2] 


mphi= 


[0, 


0] 


z= 


2 


7.7. 


ls-2 IS 


J= 





1 


43189 


nu= 


[1,2] 


spin= 


[-1/2,-1/2] 


mphi= 


[0, 


0] 


z= 


2 


7.7. 


is2s as 


J= 


1 


1 


08996 


nu= 


[1,2] 


spin= 


[+1/2,-1/2] 


mphi= 


[0, 


0] 


z= 


2 


7.7. 


ls2s IS 


J= 





1 


07930 


nu= 


[1,1] 


spin= 


[-1/2,-1/2] 


mphi= 


[0, 


-1] 


z= 


2 


7.7. 


ls2p 3P0 


J= 


2 


1 


06856 


nu= 


[1,3] 


spin= 


[-1/2,-1/2] 


mphi= 


[0, 


0] 


z= 


2 


7.7. 


ls2p 3P0 


J= 


1 


1 


06857 


nu= 


[1,1] 


spin= 


[+1/2, +1/2] 


mphi= 


[0, 


-1] 


z= 


2 


7.7. 


ls2p 3P0 


J= 





1 


06856 


nu= 


[1,1] 


spin= 


[-1/2, +1/2] 


mphi= 


[0, 


-1] 


z= 


2 


7.7. 


ls2p IPO 


J= 


1 


1 


06593 


nu= 


[1,3] 


spin= 


[-1/2, +1/2] 


mphi= 


[0, 


0] 


z= 


2 


7.7. 


ls2p IPO 


J= 


1 


1 


06594 



7. Examples for Lithium I 
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y.y. ls~2 2s 2S J=l/2 

y.y. ls~2 2p 2P0 J=l/2 

y.y. ls~2 2p 2P0 J=3/2 

y.y. ls-2 3s 2S J=l/2 



1.65595 
1.64097 
1.64097 
1.62883 



nu=[l,l,2]; spin=[-l/2, + l/2,-l/2] ; mphi=[0,0, 0]; Z=3; 
y. nu=[l,l,3] ; spin=[-l/2, + l/2,-l/2] ; mphi=[0,0, 0]; Z=3; 
7. nu=[l,l,l] ; spin=[-l/2, + l/2,-l/2] ; mphi=[0,0,-l] ; Z=3; 
y. nu=[l,l,4] ; spin=[-l/2, + l/2,-l/2] ; mphi=[0,0, 0]; Z=3; 
M=7; N=31; MAXITER=100; beta = 0; 
7. w coordinate, ranging from -1 to 1 

[Dw,w] =cheb(M); D2w=Dw~2; Hw=-diag(l-w. *w) *D2w+diag(2 . *w) *Dw; 
7. add zero weight to last term so we always use the same-sized vectors 

w_w=[inv(Dw(l:M,l:M)) (1, :) ,0] ; 
y, r coordinate, ranging from 1 to -1, rp from 1 to 0, 
y, rs from infinity to (see below) 

[D,r] = cheb(N); rp=0.5*(r+l); D = 2*D; D2 = V)~2; 

hh=diag(l-rp.*rp) ; D2=hh* (hh*D2+diag(-2*rp) *D) ; 

wr=inv(D(l:N,l:N)) (1, :) ; 
y, drop infinite weight at infinity (i.e. 1) since w.f. goes to zero quickly 
y, now wr is the same dimension as D2 etc . 

wr=wr(2:N) ./(l-rp(2:N) .~2) ' ; 
y, now multiply by hh so we didn't have zero for inversion above 

D=hh*D; 

y. Radial Laplacian minus r=infinity (for potential calcuations) 

D2inf = D2(2:N+1,2:N+1) ; 
7. Construct Laplacian minus r=0 and r=infinity 

7. (set diricelet b.c. in radial direction for eigenvalue equations) 
D2 = D2(2:N,2:N) ; D= D(2:N,2:N); rp=rp(2:N); 

zoom=l/(l/116+sqrt (beta) /41) *0 . 5 ; rs=zoom*atanh(rp) ; wr=zoom*wr; 
R = diag(l./rs); R2 = diag(l . /(rs . *rs) ) ; 
y, Hamilton for radial coordinate 
Hr=-1/ (zoom*zoom) *D2-2*R ; 

[rr,ww] = meshgrid(rs,w) ; rr = rr(:); ww = ww(:); DIMEN=numel (rr) ; 
rperp2=rr . *rr . * (1-ww. *ww) ; RR=diag(l . /rr) ; we=2*pi*kron(wr ,w_w) ; 
y. Full single electron Hamiltonian 

H=kron (Hr , eye (M+1) ) +kron (R2 , Hw) +diag (beta*beta*rperp2) ; 
mmax=0 ; 

for elcnt=l :numel (mphi) -1 

for e2cnt=elcnt+l :numel(mphi) 
if (spin(elcnt)==spin(e2cnt) ) 

dum=abs(mphi (elcnt)-mphi (e2cnt) ) ; 
if (dum>mmax) mmax=dum; end 
end 
end 
end 

y. Construct Laplacians for electron potential calculation 
Linv=zeros(DIMEN,DIMEN,mmax+l) ; 
for i=0:mmax 

y, set the denominator equal to one, if it goes to zero (don't want an infinity) 
y, we will drop rperp2=0 later 

L=kron(l/ (zoom*zoom) *D2 , eye (M+1) ) -kron(R2 ,Hw) -diag(i*i . / (rperp2+ (rperp2==0) ) ) ; 
y, set the BC for rperp2=0 but add in identity so can be inverted 
y, an infinity in the last step would have given an NaN here 

if (i~=0) L=diag(rperp2==0)+diag(rperp2"'=0)*L*diag(rperp2"'=0) ; end 
y, add in the origin as a single point 

L=[ [L, kron(D2inf (1:N-1,N) ,ones(M+l,l))] ; kron(D2inf (N, 1 :N-1) , ones (1 , M+1) ) , D2inf(N,N)]; 
y, invert Laplacian and set boundary condition at origin 

Linvhold= (eye (DIMEN+1) - [zeros (DIMEN+1 , DIMEN) , ones (DIMEN+1 , 1) ] ) *inv (L) ; 
y. divide by R before and after and remove origin! 

7. don't need to set B.C. again because the solutions for mphi''=0 are already zero in the 
y. right places, and the e-val equation sets them to zero when needed 
Linv(: , : , i+l)=-4*pi*RR*Linvhold(l :DIMEN , 1 :DIMEN) *RR; 
end 

u=zeros (DIMEN, numel(nu) ) ; 
direct_old=zeros (DIMEN, DIMEN, numeKnu)) ; 
exchange_old=zeros (DIMEN, DIMEN, numel(nu) ,mmax+l) ; 
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direct_new=zeros (DIMEN , DIMEN , numel (nu) ) ; 
excliange_new=zeros(DIMEN,DIMEN,nuinel(nu) ,mmax+l) ; 
interact=zeros (DIMEN, DIMEN) ; 
°/o begin iterations 

etotold=le33; etot=l; j=l; 

while (abs((etot-etotold)/etot)>le-5 & j<MAXITER) 
7, calculate each electron wf 
etotold=etot ; 
etot=0; 

for elcnt=l : numel (nu) 

interact=zeros (DIMEN, DIMEN) ; 
for e2cnt=l : numel (nu) 
if (e2cnt"'=elcnt) 

if (spin(elcnt)==spin(e2cnt)) 

interact=interact+direct_old( : , : ,e2cnt) . . . 

-exchange_old( : , : , e2cnt , abs (mphi ( e 1 cnt ) -mphi (e2cnt) )+l) ; 

else 

interact=interact+direct_old( : , : ,e2cnt) ; 
end 
end 
end 

if (mphi (el cnt) ==0) 

[V,Lam] = eig(H+2/Z*interact) ; 
else 

Htmp=H+2/Z*interact+diag(mphi(elcnt)*mphi(elcnt) ./rperp2) ; 
'/. apply boundary condition along z-axis, remove columns and rows with rperp2==0 

[V,Lam] = eig(reshape(Htmp(find(rperp2*rperp2'~=0)) , (N-1)*(M-1) , (N-1)*(M-1))) ; 
end 

Lam = diag(Lam) ; [Lam,ii] = sort (real (Lam) ) ; 
igood = find((V(l, :) .*V(1, :)) '<(M*N)"(-2)*le-4) ; 
V = V(:,ii); uu=V( : , igood(nu(elcnt) ) ) ; 
y, add in boundary condtion for angular direction (along z-axis) , if needed 
if (mphi (el cnt) "=0) 

uu= [zeros (1 ,N-1) ; reshape (uu,M-l ,N-1) ;zeros (1 ,N-1)] ( : ) ; 
end 

1 (elcnt)=Lam(igood(nu(elcnt) ) )+2*beta* (mphi (elcnt)+spin(elcnt) -0 . 5) ; 

u2=uu.*uu; norm=we*u2; u2=u2/norm; uu=uu/sqrt (norm) ; u( : ,elcnt)=uu; normalize 

eeen(elcnt)=2/Z*we* (uu. * (interact *uu) ) ; etot=etot+l(elcnt)-eeen(elcnt) /2 ; 

printf (' J=y.d el=y.d E-val : y.l8.12f EE: y.l8.12f Etot : y.l8 . 12f \n' , j ,elcnt ,l(elcnt) , eeen(elcnt) ,etot) ; 
direct_new( : , : ,elcnt)=diag(Linv( : , : , l)*u2) ; 

for i=l:mmax+l exchange_new( : , : , elcnt , i)=diag(uu) *Linv( : , : , i) *diag(uu) ; end 
end 

j=j+l; direct_old=direct_new; exchange_old=exchange_new; 

end 

if (j==MAXITER) 

printf ('Did not converge\n' ) ; 
end 



